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Certain features of the method of characteristics are of considerable interest in relation with 
Vlasov simulation [H. Abbasi et al, Phys. Rev. E 84, 036702 (2011)]. A Vlasov simulation scheme 
of this kind can be recurrence free providing initial phase points in velocity space are set randomly. 
Naturally, less filtering of fine-structures (arising from grid spacing) is possible as there is now a 
smaller scale than the grid spacing that is average distance between two phase points. Its interpola¬ 
tion scheme is very simple in form and carried out with less operations. In our previous report, the 
simplest model (immobile ions) was considered to merely demonstrate the important features. Now, 
a hybrid model is introduced that solves the coupled Vlasov-Fluid-Poisson system self-consistently. 

A possible application of the code is the study of ion-acoustic (IA) soliton attributes. To this end, 
a collisionless plasma with hot electrons and cold positive ions is considered. For electrons, the 
collisionless Vlasov equation is solved by following collisionless phase point trajectories in phase 
space while ions obey the fluid equations. The periodic boundary conditions are assumed. Both, 
the characteristic equations of the Vlasov equation and the fluid equations are solved using the 
Leapfrog-Trapezoidal method. However, to obtain the first half-time step of the Leapfrog, the Euler- 
Trapezoidal scheme, is employed. The presented scheme conveniently couples the two well-known 
grids in the Leapfrog method. The first test of the model is an stationary IA soliton. Trapping 
of electrons is considered and the associated phase space hole is shown. Then as a non-stationary 
test, the IA soliton generation from a localized initial profile is examined. Conservation laws are 
the other benchmark tests. 

PACS numbers: 47.11.-j, 05.10.-a, 52.30.-q, 52.65.-y 

I. INTRODUCTION 

In some plasma devices, such as the Q-macliine or 
plasma discharge, ions temperature is at least one order 
of magnitude less than electrons temperature. In such 
plasma devices, Therefore, thermal effects, associated 
with the ions, are negligible. Accordingly, for the study 
of phenomena involving kinetic effects (such as electron 
trapping), one deals with solving the Vlasov equation for 
the electrons together with the fluid equations (the con¬ 
tinuity and momentum) for the ions. An example of this 


Numerical simulations of the Vlasov equation are of 
fundamental importance for the study of many nonlinear 
processes in kinetic plasmas. The knowledge of the tem¬ 
poral evolution of the distribution function has long been 
a desire of plasma physicists as well as many involved in 
many-body physics researches. Particle trapping is an 
example where the temporal evolution of particle distri¬ 
bution function has to be considered. Many researchers 
have done great efforts with some success in the numeri¬ 
cal integration of the Vlasov equation (see the Refs. @- 
(l3| and references therein). 

In the present work, an unbounded collisionless plasma 
composed of the cold positively charged ions and hot 


* Electronic address: abbasi@aut.ac.ir 


situation is the generation of ion-acoustic (IA) soliton due 
to the nonlinear decay of a localized perturbation M- 


electrons is considered under the electrostatic approxi¬ 
mation. The Vlasov equation is solved for the electrons. 
We directly follow the characteristics along which the 
distribution function of electrons, is constant (T?). We 
choose some representative phase points, in the phase 
space, that are initialized by an initial distribution func¬ 
tion. The phase points following the characteristics are 
advanced in time by a predictor-corrector method. In¬ 
terpolation is performed between the phase points and 
a fixed grid in the phase space to obtain the distribu¬ 
tion function on the grid. From the latter, all the de¬ 
sirable quantities, such as the electron charge density, is 
obtained and used in the Poisson equation. Since the ions 
are assumed to be cold, their dynamics are governed by 
the fluid equations. Matching of the two different simula¬ 
tion schemes, needs great attention and is the main mo¬ 
tivation behind this paper. First test of the hybrid code 
is about propagation of a stationary IA soliton. Elec¬ 
tron trapping has been considered in the test problem 
as the result of their nonlinear resonant interaction with 
the IA soliton. As a non-stationary experiment, soliton 
generation from an initial Gaussian profile is considered. 
The conservation of total energy and entropy would be 
the other benchmarks. In order to avoid the error, as¬ 
sociated with the periodic boundary condition, another 
version of the code, based on moving grid, is introduced. 

The paper is organized as follows. Section II deals 
with the basic equations. Section III is devoted to the 
details of hybrid algorithm. Section IV is briefly devoted 
to first the calculation of a stationary IA soliton as a 
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test problem and then introducing a non-stationary test 
problem, i.e. disintegration of a Gaussian profile into IA 
solitons. The results of the code performance is presented 
in Sec. V. The paper terminates in section VI by a brief 
conclusion. 


II. BASIC EQUATIONS 

Let us consider, the one-dimensional electrostatic sys¬ 
tem, governing on the plasma dynamics with character¬ 
istic frequency close to the plasma frequency of ions. 
Therefore, the ion dynamics is of great importance. How¬ 
ever, we assume the thermal effects associated with the 
ions is negligible. Thus, the fluid equations are conve¬ 
nient for the ionic part of the dynamics. The electrons 
are treated kinetically. That is, the Vlasov equation gov¬ 
erns the electronic part. Poisson equation is the closure. 
In summary, we have the following set of equations: 
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where f e is the electron distribution function, v e is the 
electron phase space velocity, <f) is the self-consistent elec¬ 
tric potential, a is the ratio of electron mass to ion one 
(= m e /mi ), n e (n,) is the electron (ion) density, Vi is the 
ion fluid velocity, and the following normalizations are 
used, 
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In the above, oj p i = (47re 2 ?ro/mi) 1 / 2 , A d = 

[T e /(47re 2 n 0 )] 1/2 , no is the equilibrium value of parti¬ 
cle densities when there is no plasma perturbation, c s = 
( T e /nn ) 1/2 , e is the magnitude of the electron charge, 
and T e is the electron temperature. 


III. THE MODEL 

As it was mentioned, the electrons dynamics is the 
kinetic part of this hybrid code and is governed by the 
Vlasov equation [Eq. JT])]. In order to solve it, we directly 
follow the characteristics along which f e is constant. 

Figure 1 exhibits a typical fixed grid with Nx x Nv 
grid points (the vertices of rectangulars) and nine phase 



FIG. 1: A typical phase space grid with Nx x Nv grid points 
and 9 phase points in a host cell. Aa; and Av e correspond to 
the respective grid spacings. 



Time 


FIG. 2: The relative error in the total energy for three differ¬ 
ent population of phase points. 


points (black circles) in a cell. Thus, from hereon, we 
introduce the subscript “p” and “g” to denote the quan¬ 
tity at the phase point and grid point positions in the 
phase space, respectively. As it was mentioned in Ref. 
[12 ]. the accuracy of the code depends directly on the 
number of phase points. Figure 2 depicts the relative 
error in the total energy for three different cases. The 
curves are the results of the IA soliton experiment which 
is explained in Sec. V. Obviously, for larger population 
of phase points, the accuracy of the code is enhanced. 
Accordingly, we put nine phase points in each cell that 
is initially set regularly along X axis and to prevent the 
recurrence effect, randomly along V e axis (Fig. 1). Each 
phase point is by definition characterized by its position 
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x p and its velocity v p , and has associated with it a dis¬ 
tribution function value f p . As the phase points follow 
their collisionless trajectories, according to the following 
characteristic equations of Eq. m, 


dx p 

ir = v ^ 

(7) 

dv p E p 

dt a ’ 

(8) 


f p remains unchanged, where E p is the electric field at the 
phase point position. As the representative phase points 
follow their characteristics, they continually exchange in¬ 
formation with the fixed background grid. Each repre¬ 
sentative phase point contributes its distribution function 
to the corners of its instantaneous host cell. In this way, 
the grid distribution function, f g , is calculated from f p 
by the ” average interpolation scheme” introduced in Ref. 
0. 


space and allocating to each of them a distribution func¬ 
tion value, we know x p , v p , and f p for the Vlasov part of 
simulation. For calculating the initial value of the electric 
potential, we proceed as follows. First, by the interpo¬ 
lation fg is computed from f p . Then, by integrating f g 
with respect to velocity on the grid, is obtained. Hav¬ 
ing n° e and n°, one can solve Poisson equation by the 
well-known fast Fourier transformation (FFT) to obtain 
4>g (and therefore E g = —dfP g /dx). For this purpose, 
periodic boundary condition is assumed. (jr g can now be 
exploited in the fluid equations [Eqs. © and (dJ)], while, 
E g should be interpolated (by a third order Lagrange 
polynomial interpolation scheme jl5| ) to the position of 
phase points to obtain E p . At this stage, we have 


E° p , 


n 4 


(9) 


A. The initial loop of the code 

Let us specify each quantity at time t = nAt by a su¬ 
perscript “n”. After setting the phase points on the phase 


Our goal is to find the above quantities at t = At 
(n = 1). To do that, the Leapfrog scheme is employed 
for Eqs. Q, (JH), (Q. and (|8]) as follows 
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where the subscript “j” means the quantity at the posi¬ 
tion Xj = j Ax. Obviously, we need x]J 2 , v ^ 2 , E ^ 2 , n'J 2 , 

v)J 2 , and fl' 2 . It turns out, however, that overall accu¬ 
racy of the Leapfrog method is a very sensitive function 
of the accuracy of the half-time step quantities. In order 


to minimize the total errors, the half-time step quan¬ 
tities are computed using a predictor-corrector (Euler- 
Trapezoidal) method. For this purpose, first, we have to 
determine all quantities by the Euler method (predictor 
part) at t = At/2, 
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Now, we have the following quantities with first order of 
accuracy with respect to At, 


rV 2 
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i/ 2 
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Then, upon interpolating f p (xp^ 2 ,v^ 2 ), f^ 2 and 
1/2 

therefore n e are obtained. With the Poisson solver 
4> l J 2 (and therefore E X J 2 ) can be calculated from n^ 2 


and nj 2 . Finally, E X J 2 should be interpolated to the 

1/2 

position of phase points to obtain E p . 


The corrector scheme might be built by integrating 
Eqs. ©, ®, 0, and © using the Trapezoidal inte¬ 
gration scheme which its accuracy is second order with 
respect to At. The results is as follows, 
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The corrector part is performed in an iterative loop to 
decrease the Euler error up to a favorite value (e). Thus, 
all improved quantities after the Trapezoidal steps [Eqs. 
CHD- dMD ] of table I are compared with their previous cor¬ 
responding values and the differences can be iteratively 
reduced. A typical relative error of this kind for the ion 
velocity versus number of iterations is sketched in Fig. 3. 


This figure illustrates that performing corrector part 
is quite worthwhile. It is clear that 10 iterations are 
often enough for the relative difference of the order of 
10~ 6 . In our case 20 iterations have been used. Since, 
the initial loop is used only once in the code, the number 
of iterations is not a matter. 


J 


The outline of the procedure is given in Table I. 
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FIG. 3: The result of employing the Euler-Trapezoidal 
scheme. The typical relative error in the ion velocity cal¬ 
culation versus number of iterations. 


TABEL I. The initial loop of the code 


Initially we have: 2 $, v%, n", ip 

1. Interpolate f P {x p ,Vp) to obtain /° and n° e . 

2. Solve Poisson equation to obtain <jp g and E g from rP e and n°. 

3. Interpolate Eg to obtain Ejj. 

4. The Euler step one: determine n]/ 2 , v ]^ [Eqs. (fill) and (fl5l) ]. 

5. The Euler step two: determine xj 2 , v^ 2 [Eqs. (fTTH) and (H7[) ]. 

6 . Interpolate f p {x^ 2 , v^ 2 ) to obtain f\^ 2 and nj 2 . 

7. Solve Poisson equation to obtain 2 and E X J 2 from nj 2 and n ,/ 2 . 

8 . Interpolate E X J 2 to obtain E p ^ 2 . 

9. Trapezoidal step one: determine the “improved” n * 1 ^ 2 , v*^ 2 [Eqs. (fTfll) and (f2CT|) ] . 

10. Trapezoidal step two: determine the “improved” x* p ^ 2 , v* p ^ 2 [Eqs. ( 1211 ) and 


11. If 
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> e, and 


*1/2 1/2 


1/2 * 1/2 1/2 * 1/2 j , a 

p = x p , Vp = v p , and go to o. 
Otherwise, pass n/ 2 , h/ 2 , a^ 2 , and v p ^ 2 to the Leapfrog loop. 


> e 


B. The main loop of the code 

So far, we have defined the initial conditions and prop¬ 
erly calculated the half-time step of the Leapfrog scheme. 
Another peculiarity of the Leapfrog scheme that has to 
be noted is related to the two uncoupled grids defined in 
the Leapfrog scheme that might cause the two grids drift 
apart [16j|. In order to avoid such a decoupling of the 
grids, we proceed as follows. 

Let us first, push the phase point velocities, c™, one 
At, 

v ; +1 = v ; - — e; +1 / 2 . ( 23 ) 

Then, push the phase point positions, aip +1 ^ 2 , one At, 
x ; +3/2 = x; +1/2 + Atv; +1 . ( 24 ) 


* n _i_3 / 2 

Now, we couple the two grids. Calculating f g by the 
interpolation of /* = f p (x p +3 ^ 2 , v™ +1 ) (from hereon, the 
asterisk superscripts denote the temporary quantities). 
The temporary quantities would be corrected in the next 
steps. Next, the electron density, n * n+3 / 2 ; is computed 
by integrating with respect to velocity, 


<” +3/ 2 = J fgdv. (25) 


Now, the ion velocity, u™, is advanced one time step, 
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Then, the ion density, 


n+1/2 

Ua 


is advanced one time step 


(coupling of the two grids again), 
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The n ” +1 is the average of n ” +1 ^ 2 and n* n+3 ^ 2 , that is, Having n" +1 , one can correct n*" +3 ^ 2 , 
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In order to continue, Poisson equation should be solved 
with n* e n+3/2 and n ™ +3 / 2 that leads to 0*™+ 3 / 2 and after 
the interpolation to Ep n+3 ^ 2 . 

Now, we push the phase point velocities, v p +1 , one 
time-step, 


,j*n+2 _ yfi+1 _ ^ J^*n+ 3/2 


(30) 


Interpolation of f p (xp +3 ^ 2 , v p +3 ^ 2 ) leads to the cor¬ 
rected fg +3 ^ 2 and therefore n" +3 + Having nH +3 ^ 2 and 
n " +3 + Poisson equation can be solved to obtain </>+ 3 ^ 2 
and Ep +3 ^ 2 (after the interpolation). 


The corrected v p +3 ^ 2 is the average of v p +1 and v* n+2 , 
that is, 


v. 


"+ 3 / 2 = i « +1 + V * p n+2 ). 


(31) 


Now, the ion velocity, i>" +1 , is advanced one time step 
(coupling of the two grids one more time), 


«” +2 ) 7 = « +1 ) 1 - 


At 
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The corrected v™ +3 ^ 2 is the average of f ” +1 and v* n+2 , The outline of the algorithm of the main loop is given 
that is, in Table II. 


+ +3/2 = 2 K" +1 + <" +2 )- 


(33) 
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TABEL II. The main loop of the code 

Initially we just need: v%, Xp +1/2 , v™, v™ +1/2 , n" +1/2 , $) +1/2 , Ep +1/2 

1. Push the phase point velocities, v™, one At, to obtain v™ +1 [Eq. (1751) ]. 

2. Push the phase point positions, Xp +1 ^ 2 , one At, to obtain Xp +3 ^ 2 [Eq. (1751) ]. 

3. Calculate f* corresponding to f* = f p (xp +3 ^ 2 ,Vp +1 ) 

4. Calculate n * n+3 / 2 [Eq. (1751) ]. 

5. Advance vf, one At, to obtain v^ +1 [Eq. (j26]) J. 

6 . Advance n™ +1 ^ 2 , one At, to obtain n *" +3 A [Eq. ([ 271 ) ]. 

7. Determine n” +1 [Eq. (f28l) ]. 

8 . Determine n™ +3 ^ 2 [Eq. 


9. Solve Poisson equation with n * n+3 / 2 and n" +3 ^ 2 to obtain ^* n+3 / 2 and Ep n+3 ^ 2 . 

10. Push the phase point velocities, v™ +1 , one At, to obtain v* n+2 [Eq. (1301) ]. 

11. Determine Vp +3 ^ 2 [Eq. (1551) ]. 

12. Interpolate f p {xp +3 ^ 2 , Vp +3 ^ 2 ) to obtain fg +3//2 . 

13. Calculate nf +3 ^ 2 ■ 


14. Solve Poisson equation with n" +3 ^ 2 and n" +3 ^ 2 to obtain (f>g +3 ^ 2 and Ep +3 ^ 2 . 

15. Advance u" +1 , one At, to obtain v* n+2 [Eq. (1571) ]. 

16. Determine v™ +3 ^ 2 [Eq. (1551) ]. 

17. Pass Vp +1 , Xp +3 ^ 2 , u” +1 , u” +3/2 , n" +3 ^ 2 , </>g +3 ^ 2 , Ep +3 ^ 2 to the next step. 


IV. TEST OF THE MODEL 


Thus, we obtain, 


In order to to test the model, we examine the hybrid 
code both in stationary and non-stationary stages, let us 
first construct the stationary IA solution of Eqs. m-m- 
The stationary stage of Eqs. ([5|) and f4]) are as follows, 


d d . 

-uo—rii + — ( riiVi) = 0 , 

d£ df 

d d d 

-U 0 — Vi + Vi—Vi = <j 

d£ dll df 


(34) 

(35) 


where £ = x — uot and uq is the soliton velocity. 

Integrating Eqs. (1551) and (1551) and taking into account 
the necessary conditions for the localized profiles as £ —> 
oo 


n e ^i —> 1, </>—>■ 0, d<j)/dE —> 0, Vi —> 0. (36) 


rii = 



Vi =u 0 - 



(37) 

(38) 


The electron density is obtained from Eq. ([7]). There¬ 
fore, the stationary solution of Eq. o has to be intro¬ 
duced. Based on the polarity of soliton when —cj) is a po¬ 
tential well, a number of electrons might be in resonance 
with it and through a nonlinear mechanism, are trapped. 
The model distribution function, containing both the free 
and trapped electrons in the Maxwellian plasma was first 
introduced by E3- It is a distorted Maxwellian that has 
a hole-like structure near the IA soliton velocity, u o, as 
follows, 


// = 


a/ (27t) exp 


a/ a/ (27t ) exp 


(y/au 0 - 

— \ (\foiUQ + \/2ee) 


V e < Uq \j2(j)/c 


v e > u 0 + yj2 (j)/a 


ft = \Ja/(2it) exp -aul - /3e^J , u 0 - \/2cj)/a < v e < u 0 + y/2<j>/a 


where // and f t are the free and trapped parts of electron 


distribution function, respectively, and 
e e = ^a ( v e - u 0 ) 2 - 4>. 


(39) 

(40) 


(41) 


















FIG. 4: The stationary IA soliton that is used as the simula¬ 
tion test. 


Then, the electron density is obtained by integrating 
the distribution functions over the corresponding velocity 
range. 

Having n e and rii, the last step for the stationary IA 
soliton is the solution of Poisson equation, 


— =n e - m. (42) 

The localized solution of Eq. (1551) is numerically ob¬ 
tained and shown in Fig. 4. Equation (1451) is a nonlinear 
boundary value problem. Discretizing Eq. (1451) in the 
configuration space, we obtain 


4>j+ 1 — 2 <f>j + 4>j~i + (Ax) 2 




(43) 


Thus, the obtained potential (f> is second order accurate 
in (= Ax). In this way, we have transformed Eq. (1451) 
to a set of Ax simultaneous nonlinear equations (there 
are Nx grid points along the x axis). Recall that the 
boundary condition was introduced through Eq. (15511 . 
The set of equations have been solved by the Newton 
iterative method [lU . Each iteration deals with a matrix 
equation (containing a tridiagonal matrix) that has been 
solved by the recurrence method jl6| . 

Now it’s time to set up the non-stationary experiment, 
i.e. disintegration of an initial Gaussian profile into an 
IA soliton train. To this end, the initial potential </> is 
defined as, 


<f> = A exp 


x — C 


2 


(44) 


where A is the amplitude, A is the half of the width 
and C is introduced to fix the position of the maximum 
of the initial Gaussian profile in the simulation box. In 


order to present a variety of test problems, the polarity 
of the potential is chosen in such a way that electrons 
feel an obstacle. Therefore, the initial electron velocity 
distribution is considered as follows, 



Accordingly, the initial electron density is 

n e = exp (</>). (46) 

Then, the initial ion density can be obtained by Poisson 
equation in the following form, 


rii = n e 


d 2 (j) 

~ch?' 


(47) 


In order to define the appropriate initial ion velocity we 
proceed as follows. Since the initial profile is supposed 
to be disintegrated into several solitons through a slowly 
varying dynamics on the time scale associated with u P i 
and for each solitons Eqs. (EU) and (1551) are necessary, an 
appropriate candidate for the initial ion velocity might be 


Vi = 


2 <j>(rii - 1) 

rii + 1 


(48) 


that is obtained after substituting uq from Eq. E3 in 
Eq. (|38]l . Although, Eq. (l48l) is held in the stationary 
state for solitonic potential, our insight about the slowly 
varying dynamics led us to deduce that if the initial con¬ 
dition fulfills Eq. (1481) . the non-stationary evolution of 
the Gaussian profile will be given rise to a soliton train. 


V. THE EXPERIMENT 

In this section, the hybrid code is examined by the test 
problems introduced in the previous section. First is the 
propagation of a stationary IA soliton. For this purpose, 
assume f3 = —0.5 and uq = 1.5. Then, f p is built by 
substituting <j> (the numerical solution of Eq. [551) ] in 
Eqs. (1551) and (1451) . Moreover, the initial ion density and 
velocity are constructed by substituting <f> in Eqs. E3 
and (l38l) . respectively. Therefore, to complete the initial 
conditions, we just need to set the phase points in the 
phase space (look at the first row of Table I). In this 
experiment the following parameters are considered, 

Lx = 50, Ax = 0.05, Av e = 0.1, 

At = 0.01, V e}ma x = 300, U e ,rom = -300. (49) 

where, “Lx” is the total length of the configuration space. 
It is clear that the total length of the velocity space is 
600. Due to computational constraints, velocity cutoffs 
as in PIC models are imposed (—300 < v e < 300) [l^ |. 
Moreover, there are 9 phase points in each cell of the 
phase space, Fig. 1. 
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FIG. 5: The propagation of the IA soliton from t = 0 to 
i = 30. 


FIG. 7: The space-temporal evolution diagram of the electric 
potential. The straight line is an indication that the soliton 
moves with constant velocity. 
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FIG. 6: The electron distribution function in the phase space. 
The electron hole, associated with the IA soliton, moves to¬ 
ward the right side. 



FIG. 8: The relative error in the total energy, a) The energy 
conservation when the first step of the Leapfrog scheme is 
obtained by the Euler method, b) The energy conservation 
when the first step of the Leapfrog scheme is obtained by 
the Euler-Trapezoidal. The periodic behavior is due to the 
periodic boundary condition. 


Figure 5 depicts the results associated with the soli¬ 
ton propagation from t = 0 to t = 30. As it is seen in 
Fig. 5, the electric potential, </>, moves in the simulation 
box while its shape remains unchanged. Since /3 is neg¬ 
ative, the electron distribution function contains a hole 
in the phase space. Fig. 6 shows the electron distribu¬ 
tion function with the mentioned hole structure. It is 
obvious that it moves along the configuration axis with 
constant velocity. Note that the velocity of the hole is 
constant since there isn’t any displacement along the ve¬ 
locity axis. A space-temporal evolution diagram of the 
electric potential, in contour form, is shown in Fig. 7. In 
the case of soliton with constant velocity, the maximum 
of the electric potential has to he on a straight line of 


the space-temporal plot that its slope is equal to soliton 
speed. The figure confirms the constancy of the soliton 
velocity. The slope of the lines of Fig. 7 has been cal¬ 
culated and is 1.499. Besides, in the figure, the periodic 
boundary along the “x" axis is realizable as the repeated 
structures. 

Conservation laws are the other tests that are used to 
demonstrate the accuracy of the simulation model. The 
most basic tests are the energy and entropy conserva¬ 
tions. The model is collisionless and therefore both the 
total energy (the energy of field and particles) and the 
entropy are constants of motion. Figure 8 exhibits the 
relative error in the total energy [i.e. (total energy ra -total 
energy°)/total energy 0 xlOO, recalling that the super¬ 
script ”n” denotes the quantities at t = nAt\. As it was 
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FIG. 9: The relative error in the entropy. Whenever, the 
soliton passes the boundary, a jump takes place. 


mentioned, the Euler scheme is not enough accurate to 
be used, as the initial half-time step, by the Leapfrog 
scheme. Figure 8a demonstrates this fact as the result of 
the hybrid code execution. As it seen, at t = 100 there 
is almost three percent error. Figure 8b shows the result 
of applying the Euler-Trapezoidal scheme (see Table I) 
as the initial loop. It is clear that the error has been 
remarkably reduced. This justifies our insistence on in¬ 
troducing the Euler-Trapezoidal scheme for the initial 
loop. Figure 9 depicts the relative error of the entropy 
(that is calculated in the same manner as was used in the 
energy case). It is clear that during the execution time, 
the relative error of the entropy is limited to maximum 
0.02%. 

With a glance in Figs. 8 and 9, a kind of periodic be¬ 
havior is distinguishable. It is an important effect that 
is related to the periodic boundary condition. There are 
different way of constructing periodic boundary condi¬ 
tion (connecting the beginning and the ending of the 
configuration space). The influence of using periodic 
boundary condition exhibits itself as an inhomogeneity 
in the configuration space. Therefore, whenever the soli¬ 
ton passes the boundary, a small part of it would be 
reflected. The reflected part is very small, however, it 
would be a problem for long-time execution. That is, 
passing the boundary causes a jump in the error (both 
in the energy and entropy). To avoid this problem, an¬ 
other version of the code has been designed in which all 
the velocities has been transferred to the measure of the 
IA velocity (Galilean Transformation). Therefore, in this 
version the soliton have to be immobile and wouldn’t pass 
the boundaries. The results is shown in Figure 10. The 
figure is another indication that the soliton velocity is ex¬ 
actly u o = 1.5 and it moves without any change. Since, 
after elapsing of 100 unit of time, there is not any con¬ 
siderable difference in comparison to the initial soliton. 



FIG. 10: The immobile IA soliton as the result of executing 
the hybrid code with moving grid. 



FIG. 11: The relative error in the total energy of the simula¬ 
tion in moving grid. The reduction of error, in comparison to 
Fig. 8, is obvious 


Moreover, it is obvious from Fig. 11, that the relative 
error in the total energy has been considerably reduced 
and the periodic behavior is not seen anymore. 

The second experiment is devoted to disintegration of 
a Gaussian profile into IA soliton train. Figure 12 shows 
the result of the evolution of ion density for A = 0.2, 
A = 20, and C = 64 at T = 900. The disintegration of 
the initial Gaussian profile leads into three IA solitons 
and a linear IA wave in the back of the initial profile. 
The dotted straight line that is used to connect the max- 
imurns is an indication of well-known fact that IA solitons 
velocity (in the absence of trapped electrons) is directly 
proportional to their amplitudes. 
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T= 900 



FIG. 12: Non-stationary experiment of the Hybrid code. The 
disintegration of an initial Gaussian profile into three IA soli- 
tons and a linear IA wave. 

VI. CONCLUSION 

A hybrid scheme was introduced for simulating the 
coupled Vlasov-fluid-Poisson system. It was designed for 


a collisionless plasma with hot electrons and cold positive 
ions. The Vlasov equation was solved for the electrons 
and the ions followed the fluid equations. The periodic 
boundary conditions were assumed. The method of so¬ 
lution of the collisionless Vlasov equation was based on 
following fixed collisionless phase point trajectories. It 
was done by solving the characteristic equations of the 
Vlasov equation. Using the average interpolation scheme 
in phase space, the electron distribution function was 
mapped to a fixed background phase space grid while 
retaining it at the phase point. Both, the characteris¬ 
tics equation and fluid equations were solved using the 
Leapfrog method. However, to obtain the first half-time 
step of the Leapfrog, the Euler-Trapezoidal scheme was 
introduced. The presented scheme conveniently coupled 
the two well-known grids in the Leapfrog method. The 
first test of the model was the propagation of an sta¬ 
tionary IA soliton. The simulation code preserved the 
stationary soliton features. Conservation laws were the 
other benchmark tests. The error in the relative entropy 
and total energy was kept to less than one percent. As 
the non-stationary test, disintegration of a Gaussian pro¬ 
file into IA solitons was considered and confirmed the 
appropriate performance of the hybrid code once more. 
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